# Campionamento regolare
fs <- 1000 # freq. di campionamento
Ts <- 1/fs # periodo di campionamentoEsercitazione 1
Analisi dell’onda pressoria arteriosa
Per la consegna dell’esercitazione si faccia riferimento a questo link. Mentre le librerie sono rimandate nell’icona Librerie.
Parametri di campionamento che verranno considerati nell’analisi:
1 Costruzione onda pressoria
Un metodo per interpolare dei punti è mediante una funzione spline. Peraltro, questa permette di costruire una funzione “smussata”, priva di punti singolari nel dominio interpolato.
Inizialmente è stata costruita una spline cubica naturale, che però introduce un nuovo punto di minimo in corrispondenza di \(t \approx 0.75\,\mathrm{s}\). Dunque, per garantire che venga seguito il profilo dei punti interpolati (i.e.: è necessario seguire il più possibile l’interpolazione lineare tra i diversi punti), si possono seguire due approcci:
infittire i punti da interpolare laddove c’è una deviazione dal comportamento desiderato, ed eventualmente attribuirgli maggior importanza;
Suggerimento ChatGPT: costruire una spline monotona shape-preserving, che consente una monotonia locale al variare dei valori del dominio e quindi permette di conservare il profilo desiderato.
Come soluzione è stato scelto l’approccio proposto dall’intelligenza artificiale.
Funzione di interpolazione.
# Punti da interpolare
tind <- c(0, 0.1, 0.2, 0.3, 0.5, 1) # istanti "chiave"
tval <- c(80, 120, 100, 105, 90, 80) # valori ai tempi chiave
# Funzione spline cubica naturale
p_spline <- splinefun(x = tind, y = tval, method = "natural")
# Segnale di onda pressoria
Ta <- max(tind) - min(tind) # tempo di acquisizioni
N <- fs*Ta # numero di campioni
sig_art <- tibble(
n = 0:(N-1),
t = n/fs,
p_art = p_spline(t)
)
p1 <- sig_art %>% ggplot(aes(x = t, y = p_art)) +
geom_line(color = "blue") +
labs(x = "Tempo (s)", y = "Pressione (mmHg)", title = "Spline cubica naturale")# Spline monotona shape-preserving
p_spline <- splinefun(x = tind, y = tval, method = "monoH.FC")
# Segnale di onda pressoria
sig_art <- tibble(
n = 0:(N-1),
t = n/fs,
p_art = p_spline(t)
)
p2 <- sig_art %>% ggplot(aes(x = t, y = p_art)) +
geom_line(color = "blue") +
labs(x = "Tempo (s)", y = "Pressione (mmHg)", title = "Spline monotona shape-preserving")
p1+p2
Simulando dieci battiti cardiaci, il segnale complessivo risulta come segue.
# Ripetizione per 10 battiti
Ta <- 10
N <- fs*Ta
sig_art <- tibble(
n = 0:(N-1),
t = n/fs,
p_art = rep(sig_art$p_art, times = 10)
)
plot_ly() %>%
add_lines(x = t, y = p_art, line = list(color = "red", width = 2.5), name = "Interp. lineare") %>%
add_lines(x = sig_art$t, y = sig_art$p_art, line = list(color = "blue", width = 1.5, dash="dashdot"), name = "Spline") %>%
layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Onda pressoria")Si mette in rilievo la presenza di una cuspide ad ogni intero di periodo \(T = 1\,\mathrm{s}\). Per sistemare il problema si potrebbero imporre delle condizioni sul valore e sulle derivate tra le ripetizioni in corrispondenza di quei punti. Ciononostante, considerando che non è immediato il codice e che il segnale analizzato sarà quello ricostruito mediante le serie di Fourier, si trascura questo aspetto.
1.1 Ricostruzione mediante serie di Fourier
Il segnale viene ricostruito mediante le serie di Fourier, in cui l’\(n\)-esimo elemento vale:
\[ p[n] = a_0 + \sum_{k = 1}^{K} a_k \cos[2\pi k f_0 \, n] + \sum_{k = 1}^{K} a_k \sin[2\pi k f_0 \, n] , \]
dove \(f_0\) è la frequenza fondamentale, \(K\) è il numero totale di componenti armoniche e \(a_k, b_k\) sono le \(k\)-esime componenti, risultato della scomposizione su basi orto-normali. Si osservi che la componente in DC non viene divisa per 2, perché è calcolata come la media del segnale, anziché il doppio della media (come da teoria della trasformata).
La \(k\)-esima componente avrà frequenza pari a \(f_k = k f_0 = \frac{k}{T_a}\), dove \(T_a = 10\,\mathrm{s}\) è il tempo di acquisizione. Perciò, in virtù del Teorema di Nyquist-Shannon, deve valere la seguente disuguaglianza, per la quale si calcola la massima componente:
\[ f_{max} = \frac{k_{max}}{T_a} \leq \frac{f_s}2 \Leftrightarrow k_{max} = \frac{f_s T_a}{2} = 5000, \] con \(f_s = 1000\,\mathrm{Hz}\) la frequenza di campionamento.
Ricostruzione segnale anche con i coseni.
a0 <- mean(sig_art$p_art)
# La componente DC è ortogonale alle basi seno e cose
# Per evitare errori numerici, quest'ultima viene sottratta al segnale
p_art_osc <- sig_art$p_art - a0
K <- 500 # massimo possibile: 5e3
f_fund <- 1/Ta
comps_cos <- map_dbl(1:K, \(k) {
c <- cos_k(sig_art$t, f_fund, k)
p_art_osc %ps% c/norm_ps(c)
}) %>% zapsmall()
comps_sin <- map_dbl(1:K, \(k) {
s <- sin_k(sig_art$t, f_fund, k)
p_art_osc %ps% s/norm_ps(s)
}) %>% zapsmall()
comps_sin[1:50] [1] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
[8] 0.000000 0.000000 9.340245 0.000000 0.000000 0.000000 0.000000
[15] 0.000000 0.000000 0.000000 0.000000 0.000000 3.153407 0.000000
[22] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
[29] 0.000000 3.164646 0.000000 0.000000 0.000000 0.000000 0.000000
[36] 0.000000 0.000000 0.000000 0.000000 1.965583 0.000000 0.000000
[43] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
[50] -0.061118Routine mediante reduce.
# Implementazione mediante reduce
sig_art <- sig_art %>% mutate(
p_art_k = a0 + reduce(1:K, \(acc, k) {
c <- cos_k(t, f_fund, k)
s <- sin_k(t, f_fund, k)
acc + comps_cos[k] * c/norm_ps(c) + comps_sin[k] * s/norm_ps(s)
}, .init = rep(0, n()))
)
plot_ly() %>%
add_lines(x = sig_art$t, y = sig_art$p_art, line = list(color = "blue"), name = "Originale") %>%
add_lines(x = sig_art$t, y = sig_art$p_art_k, line = list(color = "red"), name = "Ricostruito") %>%
layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Ricostruzione del segnale")2 Filtraggio del segnale
Al segnale viene addizionato del rumore bianco gaussiano.
sig_art <- sig_art %>% mutate(
p_art_n = p_art_k + rnorm(length(n), 0, 2)
)
pp <- plot_ly() %>%
add_lines(x = sig_art$t, y = sig_art$p_art_n, line = list(color = "blue"), name = "Seg. con rumore") %>%
add_lines(x = sig_art$t, y = sig_art$p_art_k, line = list(color = "red"), name = "Seg. nominale") %>%
layout(xaxis = list(title = "Tempo (s)"), yaxis = list(title = "Pressione (mmHg)"), title = "Segnale ricostruito con rumore")
ppLa scelta della frequenza di taglio per il filtraggio può essere determinata sulla base dello spettro del segnale: trascurando la componente DC, che è energeticamente predominante ma non informativa sulla dinamica, è possibile valutare il taglio in funzione del contributo energetico associato a ciascuna componente in frequenza.
Tale approccio è giustificato dal Teorema di Parseval, il quale afferma che l’energia del segnale nel dominio del tempo è equivalente alla somma dei quadrati dei moduli dei coefficienti della trasformata di Fourier. In forma discreta, coerente con l’utilizzo della FFT, vale infatti: \[ \int_{0}^T s^2(t) \, dt = \sum_{k=1}^N |Y_k|^2. \]
Nello specifico il contributo relativo è calcolato come \(e_k = \frac{|Y_k|^2}{\sum_k |Y_k|^2}\).
p.fft <- sig_art %>% select(-p_art, -p_art_k) %>%
rename(p = p_art_n) %>%
mutate(
f = n/Ta,
fft = fft(p), # FFT concessa perché numero intero di periodi (altrimenti applicare windowing)
mod = (2-(f==0)) * Mod(fft)/length(n), # riscalatura per numero di campioni e spettro monolatero
phase = Arg(fft)/pi*180
) %>% slice_head(n = floor(N/2)) # spettro monolatero
p.fft %>% plot_ly(x = ~f) %>%
add_segments(xend = ~f, y = 0, yend = ~mod, line = list(color = "#F8766D", width = 1)) %>%
add_markers(y = ~mod, marker = list( color = "#F8766D")) %>%
layout(
xaxis = list(title = "Frequenza (Hz)"),
yaxis = list(title = "Ampiezza (mmHg)"),
title = "Spettro monolatero"
)p1 <- p.fft %>% filter(f <= 10) %>% # ci si concentra solo sulle frequenze influenti
ggplot(aes(x=f, y=mod)) +
geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#F8766D") +
geom_point(colour = "#F8766D") +
labs(x = "Frequenza (Hz)", y = "Ampiezza (mmHg)")
p2 <- p.fft %>% filter(f<=10) %>% slice(-1) %>% # si esclude la componente DC
mutate(en_rel = mod^2 / sum(mod^2)) %>% # calcola il contributo energetico relativo
ggplot(aes(x=f, y=en_rel)) +
geom_spoke(aes(angle = pi/2, y=0, radius = en_rel), linewidth = 0.2, colour = "#56B4E9") +
geom_point(colour = "#56B4E9") +
labs(x = "Frequenza (Hz)", y = "Contributo energetico")
p1/p2
Frequenza di taglio opportuna.
Osservando i grafici soprastanti, si evidenzia che la componente a frequenza \(f = 1 \,\mathrm{Hz}\) contribuisce per il 69\(\,\%\), mentre le altre sono al di sotto del \(10\%\). Pertanto, si può provare a scegliere una frequenza di taglio di poco superiore a quella precedentemente citata, ossia pari a \(f_c = 2\,\mathrm{Hz}\).